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Theoretical characterization of excitation energy transfer in chlorosome light-harvesting antennae 
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Chlorosomes are the largest and most efficient natural light-harvesting antenna systems. They contain thou- 
sands of pigment molecules - bacteriochlorophylls (BChls)- that are organized into supramolecular aggregates 
and form a very efficient network for excitonic energy migration. Here, we present a theoretical study of exci- 
tation energy transfer (EET) in the chlorosome based on experimental evidence of the molecular assembly. Our 
model for the exciton dynamics throughout the antenna combines a stochastic time propagation of the excitonic 
wave function with molecular dynamics simulations of supramolecular structure, and electronic structure cal- 
culations of the excited states. The simulation results reveal a detailed picture of the EET in the chlorosome. 
Coherent energy transfer is significant only for the first 50 f s after the initial excitation, and the wavelike mo- 
tion of the exciton is completely damped at 100 fs. Characteristic time constants of incoherent energy transfer, 
subsequently, vary from 1 ps to several tens of ps. We assign the time scales of the EET to specific physical 
processes by comparing our results with the data obtained from time-resolved spectroscopy experiments. 
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I. INTRODUCTION 

Chlorosomes of green photo synthetic bacteria are the 
largest light-harvesting antenna systems found in nature |[ll- 
0]. Depending on the bacteria species they consist of 
bacteriochlorophyll (BChl) c, d, or e molecules that form 
supramolecular assemblies without protein scaflTolding. This 
finding motivates a strong interest in the design of artificial 
light-harvestingsystems that uses aggregates of pigment or 
dye molecules i^^ . 

In green sulfur bacteria (GSB) - the model system that we 
study - light energy absorbed by the chlorosome is transferred 
in the form of Frenkel excitons via a baseplate [|1Q|] to an ar- 
ray of Fenna-Matthews-Olson (FMO) complexes llllll , which 
in turn transfers the excitation to reaction centers, where exci- 
ton dissociation occurs. Unlike the chlorosome, the baseplate, 
FMOs, and reaction centers all contain BChls a embedded in 
a protein environment. The impressive energy conversion eflft- 
ciency of the chlorosome is crucial for the survival of bacteria 
especially at low light intensities. It has been observed that 
GSB can live in environments wher e lig ht intensity can be as 
low as a few yumol photons m"^ s"^ |[l2i[l3l] . 

The EET in chlorosomes have been studied using var- 
ious time-resolved spectroscopy techniques [|14l - l23|] . For 
later convenience, we briefly review the excitation energy 
transfer (EET) time scales measured for chlorosomes from 
Chlorobaculum tepidum (C tepidum), a species of GSB. The 
chlorosome in this species consists of BChl c molecules. One- 
color transient absorption (TA) measurement in the wave- 
length region of 750-790 nm has provided four major time 
constants of 200-300 fs, 1.7-1.8 ps, 5.4-5.9 ps, and 30-40 
ps IIITI] . The authors have assigned the first time to relax- 
ation from higher to lower exciton states, the second and third 
times to the EET within BChl c aggregates in diflTerent energy- 



transfer steps, and the fourth time component to the EET 
from BChl c aggregates to the baseplate. Qualitatively sim- 
ilar results have been obtained by two-pulse photon echo and 
one-color TA experiments [18] . A recent two-color TA study 
with 685 nm pump and 758 nm probe has also yielded four 
time measurements of 120 fs, 1.1-1.2 ps, 12-14 ps, and 46-52 
ps 1 22], while the authors have regarded all time components 
as the EET within BChl c aggregates. More recently, two- 
dimensional (2-D) electronic spectra have been measured to 
explore exciton dynamics during the first 200 fs after excita- 
tion 1I23I1 . The authors have shown that the amplitude of the 
positive peak decays with an eflTective time constant of 40 fs 
and have explained changes in the 2-D spectra by proposing 
exciton diflfusion on a sub- 100 fs time scale. The experiments 
above have suggested that there are several EET processes in 
the chlorosomes from sub- 100 fs to several tens of ps, while 
interpretation of them is still controversial. 
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Electron microscopy (EM) [|24l - l26|] supported by theoret- 
ical studies 1 27] have suggested that the BChls are orga- 
nized into rod-like elements - rolls, while sheet- like aggre- 
gates 1251 ^M may coexist. This uncertainty in the structural 
characterization of chlorosomes stems from the large disor- 
der in supramolecular aggregates that compose the chloro- 
some. Recent understanding of biosynthesis of BChls in C. 
tepidum ll29n has led to a well-characterized bchT?, bchQ, 
and bchf/ three-point mutant [30]. This mutant synthesizes 
chlorosomes that contain more than 95% BChl d and form 
regularly-packed and tubular- shaped elements [1261] . Ganapa- 
thy et al. [|3l|] have proposed an atomistic structure of the roll 
from the three-point mutant by using solid-state NMR and 
cryo-EM images: Syn-anti monomer stacks are basic build- 
ing blocks, and they are organized into coaxial cylinders to 
form the tubular- shaped elements. Linear dichroism spec- 
tra of individual chlorosomes [[321] have confirmed that the 
chlorosomes from the three-point mutant are less disordered 
than those from the wild type, which is compatible with the 
highly- symmetric structure proposed by Ganapathy et al. [l3ll] 
The structure of the roll is illustrated in Figure 1 . The syn- 
anti monomer stacks (Figure lb) run perpendicular to the roll 




FIG. 1: The structural model of chlorosomes from the three-point 
mutant tSlH . (a) The schematic model of the molecular assembly, 
(b) The atomistic structure of the syn-anti monomer stacks. 



and the Qy excited state, the excitation energy of which de- 
pends on the nuclear configuration. In the single exciton ap- 
proximation, the Hamiltonian for the aggregate of BChls can 
be written as follows: 



+ r(P) + VgCR), 



(1) 



where \m) denotes the state where the excitation is localized 
at m-th BChl (site) and all other sites are in the ground states. 
6^(R) represents an excitation energy of the m-th site in the 
nuclear configuration R, and VmwCR) is a excitonic coupling 
constant between m-th and ^-th BChls. Here r(P) and VgCR) 
are the kinetic energy and the ground- state potential energy 
for the nuclear coordinates, respectively. As the next step, we 
decompose the total Hamiltonian without any assumption of 
the functional form into a system and bath contributions as: 



H = Hs -\- HsB + fi, 



(2) 



symmetry axis; ^-BChl pairs form a ring with /t-fold rota- 
tional symmetry. These rings are stacked in-parallel so that 
they constitute a helical cylinder (Figure la). Several concen- 
tric cylinders are further combined with a spacing of the order 
of2.1nm|:31]. 

In our recent letter, psl] we have shown that in the chloro- 
some based on the structure proposed by Ganapathy et al. Il3lll 
the memory eff'ects associated with the environment fluctua- 
tions - dynamic disorder - assist the exciton difl'usion for a 
broad range of static structural imperfections. In this paper, 
we present a full report of our work [|33ll and compare our 
obtained time scales to those of the experimental results men- 
tioned earlier. The structures of the roll without the baseplate 
were prepared to model exciton dynamics in the chlorosome. 
We characterize the EET by combining all-atom molecular 
dynamics (MD) simulations, time-dependent density func- 
tional theory (TDDFT) excited- state calculations, and open 
quantum system approaches. Our simulations allow for a uni- 
fied description of the coherent and incoherent energy transfer 
and depict the overall time evolution of EET processes in the 
chlorosome. 

The reminder of this paper is organized as follows: Sec- 
tion II presents the model of exciton dynamics in the chloro- 
some. Section III describes details of the site energy and ex- 
citonic coupling calculations. Section IV A shows autocorre- 
lation functions and spectral densities from MD/TDDFT cal- 
culations. Section IV B presents the optical properties of the 
roll. The EET in the chlorosome is discussed in Section IV 
C to F. Section IV G summarizes EET time scales in compar- 
ison with the experimental data. We conclude the study by 
summarizing our results in Section V. 



II. THEORY 

We first introduce the excitonic Hamiltonian for the chloro- 
some. Each BChl is treated as a two-level system of ground 



^5 = ^ <^m) |m)<m| + ^ Vm«|^)<^|, (3) 



//55 = ^A6^(R)|m)<m|, (4) 



Hb = r(P) + Vg(R\ (5) 

where A6;„(R) = em(R) - {^m}- The brackets () denote the 
average over the nuclear configurations, and the Condon ap- 
proximation was used in order to omit the dependence of Vmn 
onR. 

The time dependence of the site energies em(R) comes from 
the dependence of their coordinates R on time t. This dynamic 
disorder arising from the nuclear motions is determined by 
the atomistic simulations. We use MD simulations to gener- 
ate R(0 and TDDFT excited- state calculations to obtain e(R). 
Thus, in contrast to many studies based on a quantum master 
equation, this approach can describe the system-bath coupling 
in complete atomistic detail. 

In addition to the dynamic disorder, there is another disor- 
der among excitation energies of the BChls which does not 
change on the time scales of the exciton dynamics. The phys- 
ical origins of this static disorder are due to the existence of 
diff'erent homologues of the BChls [30, M^^5J, structural im- 
perfections of supramolecular arrangements, and variations 
among individual chlorosomes. Experimentally, the disorder 
among individual chlorosomes is estimated from low temper- 
ature spectral hole burni ng [UVl 1361] and single molecule flu- 
orescence spectroscopy ||32l [37l - l39l] . It is difficult to predict 
this static disorder from MD simulations. Thus we account 
for the static disorder by introducing random shifts in the av- 
erage site energy from a Gaussian distribution f(em - {^m)) = 
1/ y27To-^Qxp(-(€m- {€m))^/2cr^\ whcrc cr is the standard 
deviation for the static disorder. The cr is assumed to be iden- 
tical for all BChls. 



If we treat the nuclear degrees of freedom classically, the 
wave function of the exciton system is described with the 
time-dependent Schrodinger equation [|4Q|] : 



in-\ilf{t))^m)\iif{t))^ 



dt 



(6) 



where H{t) = Hs -\- HsbO)- In this approach, short MD trajec- 
tories are sampled from a full MD trajectory, and the excitonic 
system is propagated under a unitary evolution for each short 
MD trajectory. Here it is assumed that the excitonic system 
does not affect the bath dynamics. The density matrix of the 
excitonic system p is obtained as an average of these unitary 
evolutions: 



1 ^ 

P(0 = -XM><^/(0|. 



(7) 



A similar approach has been previously applied to small light- 
harvesting systems such as the FMO complex ll^llfi^ and 
light-harvesting II complexes of purple bacteria MM nM • 

For a large number of BChls present in the chlorosome, 
it is computationally unfeasible to calculate site energies of 
all pigment molecules along an MD trajectory. Here, we 
take an approximate approach to obtain the site energy fluc- 
tuations. Following our previous work on the FMO com- 
plex 1I41I1 . the site energy fluctuations are obtained by solving 
a set of Langevin equations with white noise: 



dt Tm 



(8) 



where Fmif) is a random force; {Fmif)) = and 
{Ffn(t)Fm(0)} = 2^A6^)^(0/t,^. a correlation function from 
eq[8]is 



Using eq [TOl the time evolution of the density matrix can be 
obtained as follows Hill]: 



|p(0 = - ^ [Hs,p(t)] + ^rm[A.p(OAl 

m 

-^A^Alp(t)-^p(t)A^All (11) 

where A^ = \m){m\. Here, a pure dephasing rate for the m-th 
site y,n is given by fsd] n^jm = 2 (a^) t^. 

We also construct a classical hopping model using the fol- 
lowing equation [51]. 

T-;?m(0 = ^ knmPn(t) - ^ kmnPm{t), (12) 

n n 

where the hopping rate from m-th to ^-th sites kmn is given 
by 01, 






d(jjA^{(jj)Fn{oj). 



(13) 



Here, pm is a population for m-th site, and Am^co) and Fn^oj) 
are normalized absorption and emission and spectra of the m- 
th and ^-th sites, respectively. With the use of eq[TOl the ab- 
sorption and emission lines can be derived as Lorentz func- 
tions iDin so that the overlap of the absorption and emission 
spectra becomes: 



In I dajAtn(oj)Fn(oj) 

kJ — OO 



(14) 



where Hojm = {^m)- 



(A6^(0A6^(0)) = (A6^)exp 



This procedure gives the same results as the Kubo-Anderson's 
stochastic model ||4^ |47|]. The relaxation time of the 
bath fluctuations, Tm, can be determined by integrating au- 
tocorrelation functions (ACF) of the site energies, r^ = 
^ dt{A€m(t)A€m(0)} / (Ae^). The input parameters to this 

description of the exciton dynamics are {Aem), [Ae^V and Tm, 
which are obtained from the MD/TDDFT calculations. Owing 
to the symmetry of the roll, we can assume that the statistical 
properties of the environmental noise are similar for all the 
syn-anti stacked monomers. Therefore, a random sampling 
of MD/TDDFT trajectories for a single dimer of syn and anti 
BChls can be used for the whole roll. 

In order to check the memory efl'ects, we compare our 
model with the Haken-Strobl-Reineker (HSR) model HsHH] 
where the site energy fluctuations do not depend on their 
memory. In the HSR model, the correlation functions of the 
site energies are approximated as 



(A6,,(0A6^(0)) = 2(A6^)t^^(0. 



(9) III. SITE ENERGY AND EXCITONIC COUPLING 

CALCULATIONS 

We set up a model of three-concentric five-stacked (3x5) 
rings (Figure 2), so that polarizations from several layers can 
be incorporated into excited-state calculations. Rings from 
diff'erent layers consist of 60, 80 or 100 BChl molecules, and 
their radii are 6.1, 8.2, and 10.2 nm, respectively. The total 
number of the BChl molecules in the structure is 1200. 

The force field parameters of the BChl d were prepared 
as follows: One syn-anti monomer unit was extracted from 
Ganapathy's original structure, and the Mulliken charges of 
the dimer were used as partial charges in the force field. The 
Mulliken charges were calculated with the Becke three param- 
eter Lee- Yang-Parr hybrid functional 0521 153|] and 6-3 IG** 
basis set. Other parameters were adopted from the General 
Amber Force Field [1541] . 

After the geometry optimization, an MD simulation was 
performed with a time step of 2 fs. The temperature was set 
at 300 K by the Langevin thermostat. The hydrogen atoms 
were constrained at their ideal positions using the SHAKE al- 
(10) gorithm. After a 1 ns equilibration run, a 20 ps production run 



(a) 



j^^^^^ 



(b) 
(c) 




FIG. 2: (a) Top view and (b) side view of a model of three-concentric 
five- stacked rings for site energy and excitonic coupling calculations, 
(c) Syn (red) and anti (blue) monomers in the three-concentric five- 
stacked rings, the site energies and the TrEsp charges of which were 
obtained from the excited- state calculations. In (a) and (b), famesyl 
tails of the BChls are not shown for clarity. 



was carried out and stored every 4 fs. The MD simulation was 
performed by a CUDA implementation of NAMD ^^ [56ll . 
Using the snapshots from the 20 ps trajectory, the Qy excita- 
tion energies for the syn and the anti monomers (Figure 2c) 
were obtained. Then, the averages and the two-times corre- 
lation functions of the site energies were calculated. TDDFT 
excited- state calculations were performed with the long-range 
corrected hybrid functional of Becke [dtIi and 6-3 IG ba- 
sis set in the presence of external point charges of all other 
molecules, as implemented in the Q-CHEM quantum chem- 
istry package 15 8i] . Results were interpolated to 2 fs time 
steps. 

In our previous paper [1331] . we employed the point-dipole 
approximation (PDA) with experimental estimates of the tran- 
sition dipole moments. This approximation is strictly applica- 
ble only if the intermolecular distance is much larger than the 
size of the molecules |[53, [60], which is not the case in the 
chlorosome where both sizes are comparable. In order to re- 
lax the PDA, we decided to apply transition charges from the 
electrostatic potentials (TrEsp) method l,6Qi] as it offers a com- 
promise between accuracy and computational times. 

In the TrEsp method 11601] . the Coulomb interaction be- 
tween transition densities is approximated as a sum of pair- 
wise interactions of atomic charges. We employ optimally- 
weighted charges |61] to obtain partial atomic charges from 
transition densities, where Mulliken charges of the transition 



densities were used as the reference. A single pair of syn-anti 
stacked monomers in the 3x5 rings were computed quantum- 
mechanically (Figure 2c), and all other molecules were treated 
as external point charges. The transition densities were cal- 
culated using multi-layer fragment molecular orbital meth- 
ods Il62n in conjunction with configuration interaction singles 
and the 6-3 IG basis set. ABINIT-MP(X) fH,!!! program 
package was used for the excitonic coupling calculations. Fi- 
nally, the TrEsp charges of the syn and anti monomers were 
used to calculate all excitonic couplings in the structure. 



IV. RESULTS AND DISCUSSION 

A. Autocorrelation functions and spectral densities 

The results of the MD/TDDFT calculations are summarized 
in Table I. The calculated average site energies are 15243.0 
and 15319.4 cm"^ for the syn and the anti monomers, respec- 
tively, and the corresponding standard deviations are 532.2 
and 556.5 cm"^. Although the syn and the anti monomers 
do not experience the exact same environment, the difference 
in the averages and standard deviations of the site energy are 
small. The ACFs and spectral densities for the syn and the anti 
monomers are shown in Figure 3. The spectral densities were 
calculated as a weighted cosine transform of the correlation 
functions m Hill: 



JmM = 



CD r 

T^ksT Jo 



dt <A6;^(0A6^(0)) cos {ojt). (15) 



In eq [15] we used a harmonic prefactor in order to negate 
the temperature dependence of the classical correlation func- 
tion. 11641] In the obtained temperature-independent spectral 
density. Figure 3b, the strong peaks around 1600 to 2000 cm"^ 
can be attributed to internal bond stretching in the porphyrin 
macrocycle of the BChl, while the low-frequency region is as- 
cribed to intermolecular interactions arising from the closed- 
packed structure. 

The reorganization energies 658.7 and 719.9 cm"^ for syn 
and anti monomers, respectively, were calculated as Er = 
C^ d(jjj((jj)/(jj. These values are larger than those calculated 
from the spectral densities with the standard prefactor [641] in 
our previous paper [1331] . Our MD/TDDFT calculations imply 
that the reorganization energies of the BChl d in the roll is 
larger than that in solvent (69 cm"^) II65|] . This large reorga- 
nization energy can be a signature of a fast energy relaxation 
within Qy exciton states in the roll. 

The chlorosome pigments should be distinguished from 
ideal J-aggregates where the J-band with the largest oscilla- 
tor strength is also the lowest exciton state (LES) IIITI 1361] . 
Psencik et al. iflTI] has shown that LESs in chlorosomes from 
C tepidum are distributed within 760-800 nm range, while 
the absorption maximum is at 750 nm. We hypothesise that 
the relatively large reorganization energy enhances the re- 
laxation from the absorption maximum to the lower exciton 
states, leading to the energy funneling to the baseplate. 

The relaxation time constants for the syn and the anti 
monomers were calculated to be 8.80 and 4.5 1 fs, respectively. 
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FIG. 4: (a) A model of ten-stacked rings for the exciton dynamics, 
(b) Initially excited monomer (red) and its neighboring monomers in 
the ten- stacked rings. In (a) and (b), famesyl tails of the BChls are 
not shown for clarity. 



FIG. 3: (a) Autocorrelation functions (ACF) of the site energies and 
(b) spectral densities for the syn (red) and the anti (blue) monomers. 
Low-frequency regions of the spectral densities are shown in the in- 
set. 



These values are in good agreement with those estimated by 
the atomistic simulation for the FMO complex lIlllS^l. How- 
ever, the bath correlation time obtained from fitting to the 
experimental absorption spectrum of the chlorosome [23] is 
several times longer than the present values. The reason for 
this discrepancy can be, for instance, due to the fact that 
our ab initio model assumes a weak electron- vibrational cou- 
pling. Moreover, our model accounts for static disorder in the 
chlorosome, which was not included in the phenomenological 
model from Ref. 23. We expect that our calculations provide 
a lower bound estimate for the bath correlations times. 



B. Spectra 

To compute absorption and circular dichroism (CD) spectra 
we run simulations of the exciton dynamics using ten-stacked 
rings (Figure 4a), where each ring is composed of 60 BChls. 
Initial conditions of the site energy fluctuations l^^emit = 0) 
were obtained from the Gaussian distribution with a variance 
(Ae^Y and the Euler-Maruyama scheme was applied to solve 
eq[8]in order to obtain A6^(0 and define H(t). Initial condi- 
tions of the wave functions were chosen as a localized state on 



an anti monomer (Figure 4), and wave functions were prop- 
agated with the time step of 1.0 fs by diagonalizing H(t) at 
each step. The quantum trajectories were averaged over 1000 
diff'erent realizations of each static disorder. 

The absorption and CD spectra were calculated from the 
exciton propagation by using the following expressions [l4ll 
I4I: 



x{{U^„(t,0))-{u:jt,0))], 



(16) 



,)> 



Icd{oj) oc 

Xoo 
J^^^""^^<(R^-RJ-Ou^Xyu,) 
m,n 

x{{UUt,0)}-{UUt,0))], (17) 

where UmnO, 0) is an (m,n) element of the time evolution op- 
erator for the H(t). 

Treating cr as a phenomenological parameter, we first es- 
timate its value by comparing the calculated spectra with the 
experiments. Specifically, we compare our result with spectra 
measured for the wild type of bacteria because optical prop- 
erties of the three-point mutant are not well established. The 
full width at half maximum (FWHM) of the absorption spectra 
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FIG. 5: (a) Full- widths at half maximum (FWHMs) of absorption 
spectra (red-solid) and the lowest excited state (LES) (blue-solid) as 
a function of a standard deviation of the static disorder. Experimental 
linewidths of the absorption spectra (red-dotted) and the LES (blue- 
dotted) were taken from Refs. 23 and 17, respectively, (b) Peak posi- 
tions of the absorption spectra (red) and the LES (blue) as a function 
of a standard deviation of the static disorder. 



were calculated as a function of a and compared with the ex- 
perimental value from Ref. 23. As both the dynamic and static 
disorder are included in the exciton dynamics, the linewidth of 
the absorption spectra computed with eq [16] reflects both ho- 
mogeneous and inhomogeneous broadening. 

In addition to the FWHM of the absorption spectra, we 
further validate our model by comparing the inhomogeneous 
broadening of the LES. The inhomogeneous broadening of 
the LES has been experimentally estimated from the FWHM 
of the fluorescence spectrum at the temperature of 4 K wh . 
To calculate the corresponding linewidth, the LESs were ob- 
tained by diagonalizing the system Hamiltonian with the static 
disorder, and the FWHM were estimated from a standard devi- 
ation of 10000 difl'erent realizations of the static disorder. The 
resulting FWHM corresponds to the inhomogeneous broaden- 
ing of the LES. 

Figure 5a shows the FWHMs of absorption spectra and 
the LES as a function of cr in comparison with the experi- 
ments [[TtIIH. We found that the both FWHMs can be repro- 
duced with cr = 300 cm"^ 



FIG. 6: (a) Absorption and (b) CD spectra obtained from exciton 
dynamics with the static disorder cr - 300 cm~^ . 



The absorption and the CD spectra with cr = 300 cm"^ are 
shown in Figure 6. The CD spectrum of chlorosomes from 
the wild type of C. tepidum, for example, has two negative 
and one positive band 1I22I1 . whereas there is one positive and 
one negative peak in another experiment |i66] . Our CD spec- 
trum agrees with the one reported by Saga et al. [661 an d the 
recently measured spectra for the three-point mutant [1671 . Our 
simulated result further supports the supramolecular structure 
proposed by Ganapathy et al. |3 1 ] 

Figure 5b shows the peak positions of the absorption max- 
imum and the LES as a function of a. The peak positions 
are overestimated compared with the experiments mh . One 
of the reasons for that is the small basis set in the TDDFT 
excited- state calculations. The energy difl'erence between the 
absorption maximum and the LES is only 23 cm"^ without 
the static disorder. The absorption maximum is insensitive to 
cr, while the LES decreases with increasing cr. The energy 
difl'erence at cr = 300 cm"^ is 457 cm"^ and in reasonable 
agreement with about 530 cm"^ estimated for the wild type 
chlorosome [17]. 

In the following, we discuss how these results are re- 
lated to experimental observations. Experimental studies have 
shown [130, [35I] that absorption spectra of chlorosomes from 
cells grown at low light intensities exhibit red shift up to 100 
cm"^. Borrego et al. [|35ll have found that the red shift is 



correlated with content of BChl c homologues more methy- 
lated at C-8^ and C-12^ positions. Chew et al. BOO have 
discussed possible effects of those modifications of BChl c. 
One possibility is that the modifications stabilize formation of 
larger BChl aggregates which can be assembled into a large 
chlorosome containing a greater number of BChl c. An- 
other possibility is that the modifications increase disorder 
in the aggregates and allow a much higher packing density 
of BChl c in the chlorosomes without causing crystalliza- 
tions. Our theoretical calculations show that the peak posi- 
tion of the absorption maximum becomes red-shifted as the 
static disorder increases. This tendency agrees with the ex- 
perimental observations. However, the computed shift can- 
not account for the whole red-shift of the absorption line ob- 
served in experiments. Our results are also consistent with the 
larger linewidths y3m and larger Stokes shifts [l35ll observed in 
chlorosomes grown at lower light intensities. It should be no- 
ticed that we treat the static disorder with a simplified model 
of random shifts in the site energies. In contrast, methylation 
of BChl molecules will have an influence on both site ener- 
gies and excitonic couplings by distorting the supramolecular 
arrangements. It remains to be seen how those modifications 
induce the disorder in the site energies and excitonic couplings 
at a molecular level. 



C. Coherent energy transfer 
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Here, we will discuss the coherent energy transfer within 
the roll. In order to show how an exciton is delocalized af- 
ter the initial excitation of a single site, we compute exciton 
dynamics using eqs [6l [71 and [8l In addition to the site pop- 
ulations, we present pairwise concurrence between initially 
excited site and other sites to characterize the role of coher- 
ence in the exciton dynamics. The pairwise concurrence mm 
is defined as 2 |(m \p{t)\ /)|, where / and m refer to the initially 
excited site and other sites. 

Figure 7 shows the population and coherent dynamics of 
selected monomers calculated with cr = 300 cm"^. The popu- 
lation of initially excited site decays quickly with a time con- 
stant of 25 fs without a oscillatory behavior, which reflects 
large exciton- vibration couplings. The coherences increase 
and decay within 40-50 fs; these values give an estimate of a 
decoherence time. The results suggest that the exciton is co- 
herently delocalized within 40-50 fs after the initial excitation 
of a localized state on a single site, and the coherent propaga- 
tion is quickly damped due to the dynamic disorder. 

In Figures 8 and 9, we show the population and coherence 
dynamics among the roll at selected times. By comparing Fig- 
ures 8 and 9, the coherences spread and are suppressed more 
rapidly than the populations, while transport directions are the 
same. These features are qualitatively similar for all examined 
values of the static disorder. The results show that the most 
of the excitation is distributed over neighboring BChls during 
the 100 fs after initial excitation, which is comparable to the 
sub- 100 fs exciton diffusion model recently proposed in Ref. 
23. 



FIG. 7: (a) Populations and (b) coherences of the selected sites neigh- 
boring to the initially excited site, which are obtained with the static 
disorder cr = 300 cm~^ The site numbers refer to corresponding 
BChls in Figure 4b. 
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FIG. 8: Populations obtained with the static disorder 300 cm~^ at 
selected times after the initial site excitation. Pink spheres represent 
Mg atoms of the BChl molecules, and the size of the Mg atom reflects 
the population value of the BChl molecule. 
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FIG. 9: Coherences obtained with the static disorder cr - 300 cm~^ at 
selected times after the initial site excitation. Pink spheres represent 
Mg atoms of the BChl molecules, and the size of the Mg atom re- 
flects the coherent value between the BChl molecule and the initially 
excited BChl. 



D. Exciton diffusion within the roll 

To quantify exciton transport on a timescale longer than 
the decoherence time we compute second moments of the site 
populations as a function of time. The second moments of the 
exciton propagation throughout the circumference (M^l) and 
the coaxial (Mzz) directions (see Figure 4a) calculated with 
(T = 300 cm"^ are shown in Figure 10. Both components of 
the second moment initially scale quadratically and then lin- 
early in time. The transition from the initial ballistic regime to 
the diffusive regime is observed at 20-30 fs, which is consis- 
tent with the time scale of suppression of quantum coherence 
discussed in the previous section. In contrast to our earlier 
study with the PDA couplings Il33n . the results show that M^l 
is comparable to Mzz- It has been suggested that the syn-anti 
monomer stacks run perpendicular to the symmetry axis in the 
chlorosome from three-point mutant [31], while they run par- 
allel to the symmetry axis in the chlorosomes from the wild 
type of bacteria ||3Tl l69|] . Our results may imply that exciton 
diffusion is insensitive to the supramolecular arrangements, 
which is relevant to the robust EET. 

The second moments calculated using the HSR and hop- 
ping model are also shown in Figure 10. The KA model gives 
the fastest spread of exciton population over the roll, leading 
to memory-assisted diffusion (MAD). As discussed in Section 
II, the difference between the KA and the HSR is due to the 
memory of the bath fluctuations, which is included in the for- 
mer model. This comparison indicates that the memory of the 
bath fluctuations enhances the diffusion. The HSR and the 
hopping models have the same slope with the difference only 
in the initial ballistic propagation up to around 30 fs. 

The mechanism by which the memory effects enhance the 
diffusion can be explained by considering site-to- site hopping 
rates. The hopping rate is determined by the overlap be- 
tween donor fluorescence and acceptor absorption spectra. If 
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FIG. 10: Second moments of the exciton populations in the (a) cir- 
cumference {Mil) and (b) coaxial {Mn) directions (see Figure 4a) 
with the static disorder a = 300 cm~^ . Three models compared are 
Kubo-Anderson (black), Haken-Strobl-Reineker (red), and classical 
hopping (blue) models. 



the correlation function of the site energy is a delta function 
(HSR model), absorption and fluorescence lineshapes become 
Lorentzian so that the Franck-Condon factor is an overlap of 
two Lorentz functions (eq[T4]). On the other hand, if the corre- 
lation function is exponential like the Kubo-Anderson model, 
lineshapes are closer to the Gaussian [ItqI] . It follows that the 
KA model gives larger overlap between the donor absorption 
and the acceptor fluorescence spectra. Therefore, the memory 
effects of the bath fluctuations give larger hopping rates be- 
tween pigments and thus lead to the larger exciton diffusion. 
Note that the memory effects is more pronounced if the bath 
correlation time is larger than the present value as suggested 
from the experimental study [1231] . 

Finally, we investigate the dependence of the exciton dif- 
fusion on the static disorder. Diffusion coefficients in the cir- 
cumference and coaxial directions were calculated as a func- 
tion of cr and are shown in Figure 1 1 . The HSR model gives 
almost the same diffusion constants as the hopping model: 
the difference is only in the initial ballistic propagation up to 
around 30 fs. The KA models gives about twice as large diffu- 
sion coefficients as those from the HSR without the static dis- 
order, while the memory effects become small with increasing 
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FIG. 1 1 : Exciton diffusion coefficients in the (a) circumference (Dll) 
and (b) coaxial directions as a function of a standard deviation of the 
static disorder. Three models compared are Kubo- Anderson (black), 
Haken-Strobl-Reineker (red), and hopping (blue) models. 



cr. The time when an exciton arrives at the opposite site of 
the roll from the initial site excitation can be used as the time 
scale of exciton equilibration over the single-layer roll, where 
the exciton populations are equal for all sites. This time can 
be calculated as (Rc7i)^/2Dii, where Re is the radius of the 
roll. For the roll of radius Re = 6.1 nm (60 BChl molecules 
per a ring) with the static disorder of 300 cm" ^ we obtain the 
equilibration time to be 8.1 ps. 



E. Inter-layer and inter-roll energy transfer 

So far we have focused on the single-layer roll, yet chloro- 
somes may consist of several multi-layer rolls as suggested 
from EM images 11261 [3lll . Here we explore the EET among 
different layers and rolls. We set up two-concentric ten- 
stacked (2x10) rings and two-roll system consist of 2x10 rings 
for analysis of the inter-layer and inter-roll EET. Rings from 
the inner and outer layers consist of 60 and 80 BChls, respec- 
tively. Their structures and initial conditions for exciton dy- 
namics are shown in Figure 12b and 13b. In order to speed 
up our quantum dynamics simulations, we have developed a 
code for large systems which can propagate wavefunctions in 
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FIG. 12: (a) Exciton populations of the inner (red) and outer layers 
(blue) in two-concentric ten- stacked rings. The structure and exciton 
populations at (a) and (b) 1 ps are also presented. The pink spheres 
represent Mg atoms of the BChl molecules, and the size of the Mg 
atom reflects the population values. In (b) and (c), atoms surrounding 
the Mg are only shown. 



parallel. The simulation conditions are same as those in the 
earlier sections, and exciton populations were obtained with 
cr = 300 cm-\ 

Figure 12a shows the populations of the inner and outer 
layers in the 2x10 rings. The single exponential fitting of the 
inner population yields the time constant of 1 . 1 ps, and double 
exponential fitting does not improve the fit. The time of the 
inter-layer EET is faster than that of the exciton equilibration 
within the layer. The time scale of exciton equilibration over 
the outer layer can be estimated as 14.4 ps by taking into ac- 
count the diameter of the outer layer. The result suggests that 
for a localized initial excitation the exciton is first distributed 
among different layers with the characteristic time of 1.1 ps 
and then equilibrated over the roll in 8.1-14.4 ps, respectively 
for the inner and outer rolls. 

Figure 13a shows the population dynamics between the two 
double rolls. Exponential fitting for the population of the left 
roll gives two time constants of 1.4 ps (26%) and 69.8 ps 
(76%). The faster time component of 1.4 ps is compared to 
that of inter-layer transfer and thus describes the energy trans- 
fer between outer layers of the left and right rolls. Because of 
the small contact between the rolls, the inter-roll energy trans- 
fer becomes slower after the exciton diffuses over the left roll. 
The inter-roll EET transfer occurs on major time scale of 69. 8 
ps, which is much larger than the exciton equilibration over 
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FIG. 13: (a) Exciton populations of the left (red) and right (blue) rolls 
in the system composed of two two-concentric ten-stacked rings, (b) 
Top view of the structure and the populations at ps (the initial con- 
dition of the exciton dynamics). The exciton populations at (c) 0.1 
and (d) 50 ps are also presented. The pink spheres represent Mg 
atoms of the BChl molecules, and the size of the Mg atom reflects 
the population values. In (b)-(d), atoms surrounding the Mg are only 
shown. 



the single roll. Note that the two-roll system is an extreme 
structure. In chlorosomes, a roll may be surrounded by other 
rolls or lamellar layers [1261 . Thus, the estimated time con- 
stants 1.4 and 69.8 ps can be considered as the lower and the 
upper bounds of the inter-layer energy transfer time scale. 



F. Anisotropy decay 

Fluorescence anisotropy decay has been used in the analy- 
sis of the lifetimes of coherent energy transfer. [ItiI - IttIi Fol- 
lowing Ref. 77, we calculate the fluorescence anisotropy de- 
cay based on the exciton dynamics. In the exciton propa- 
gation model we have added the ground state and the ofl'- 
diagonal system-radiation Hamiltonian, Hij(t) = E(t) • jiij to 
eq O Here E(0 is an electric field for a stimulating pulse, 
and jdmn is a dipole moment operator. An (ij) component of 
polarization tensor atj is defined as the i-th component of the 
dipole moment after applying an electric field polarized in the 
j direction. The anisotropy decay r(t) can be calculated as 
r(t) = 3y^/(45a^ -\- 4y^), where a and y are isotropic and 
anisotropic tensor invariants of polarization tensor, respec- 
tively uM- An initial condition was chosen as the ground 
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FIG. 14: The polarization anisotropy decay of the three-staked rings 
with the static disorder cr = 300 cm~^ 



state p = |G)<G|, and an 8 fs Gaussian-enveloped electric field 
of frequency 14890 cm"^ was applied as the stimulating pulse. 
We use the three-stacked rings to calculate the anisotropy de- 
cay. 

Figure 12 shows the polarization anisotropy decay for 
the three-stacked rings structure. The initial value of the 
anisotropy is close to 0.8, and the residual value is about 
0.18. The anisotropy shows the oscillatory behavior within 
50 fs after the stimulating pulse, reflecting wavelike motion 
of the exciton. Single exponential fitting to the anisotropy 
profile yields 30 fs as the depolarization time. The depo- 
larization time can be compared to the decoherence time 
estimated from the coherences in site basis in Section IV 
C. Note that the ultrafast depolarization time of sub- 100 fs 
has not been resolved by earlier time-dependent polarization 
anisotropy experiments flTsl [161 S 12211 • For example, the 
one-color TA anisotropy in the wavelength region of 720-790 
nm has given a single decay time of 1.7-3.9 ps |16]. How- 
ever, the ultrafast decay of the anisotropy has been observed 
by a recent pump-probe measurement of the chlorosome from 
Chlorobium Phaeobacterioides (Jeongho Kim personal com- 
munication), which is consistent with the present result. 



G. Time scales 

In this final section, we summarize the various time scales 
of EET obtained from our simulations. Table II shows 
the theoretical time scales in comparison with experimental 
time scales measured for chlorosomes from wild type of C. 
tepidum. It should be noticed that the EET time scales may 
be aff'ected by the presence of the baseplate acting as an ex- 
citon drain, which lies outside the scope of this paper. Ex- 
citon dynamics in the whole light-harvesting system, which 
is composed of the chlorosome, baseplate, and FMO com- 
plexes, is being studied in our group, which will be submitted 
elsewhere. 

The experimental studies (itI [isl |22l have suggested that 
there are four major time constants in the EET processes. 
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Here, we compare our results with the time components of 
120 fs, 1.1-1.2 ps, 12-14 ps, and 46-52 ps which were ob- 
tained with 685 nm pump and 758 nm probe 11221] . which 
should have a minimal contribution from the baseplate BChl 
a. Our simulations show the coherent propagation of exci- 
tons for the initial 50 fs after excitation. The lifetime of this 
coherent energy transfer can be estimated from the decoher- 
ence time and the depolarization time. The first component 
of 120 fs can be associated with the coherent EET which de- 
scribes wavelike motions of the exciton. The depolarization 
times of 30 fs is smaller than the first component of the TA 
experiments. A possible reason for this is that in terms of 
the Bloch equation interpretation, the time components of the 
TA experiments will describe the population decay (Ti) time, 
while the depolarization time includes a pure dephasing time 
(T2) as well ItiI ItsIi. The incoherent EET subsequently takes 
places after the ultrafast process. We have obtained 1.1 ps for 
the inter-layer EET within a single roll and 1 .4 ps for the en- 
ergy transfer between neighboring layers of the double rolls. 
The time scale of 1.1-1.4 ps is in good agreement with the 
second time component of 1.1-1.2 ps; the second component 
describes the inter-layer energy transfer. The 8.1-14.1 ps time 
is comparable to the third component of 12-14 ps, suggest- 
ing that the third component reflects the exciton equilibration 
within single roll. The time constant of 69.8 ps corresponds 
qualitatively to the 46-52 ps. The forth time component is 
tentatively assigned as the exciton equilibration over difl'erent 
rolls, while this would not be the case if a roll were surrounded 
by lamellar layers. 

Based on the interpretations above, we speculate on the cor- 
responding time scales of the EET in the chlorosomes based 
on the structure from the three-point mutant. The first com- 
ponent describes the lifetime of the coherent energy transfer 
which largely depends on the exciton- vibration coupling. Re- 
organization energies will be similar for the wild type and 
the three-point mutant as long as the local packing structure 
is same. Therefore, the first component reflecting the coher- 
ent energy transfer may be similar for the wild type and the 
three-point mutant. We also expect that the second compo- 
nent describing the inter-layer EET is similar for the wild type 
and the three-point mutant, because the same 2.1 nm spacing 
was observed for both the wild type and the three-point mu- 
tant ll^lsil] . The third component describing the energy equi- 
libration over the roll is sensitive to a radius of the roll. The 
EM study 1I26I1 finds that the diameter of the roll of the three- 



point mutant is larger than that of the wild type. According to 
our interpretation and the experimental observation, the third 
time component of the three-point mutant will be larger than 
that of the wild type. This could be an experimental verifica- 
tion of our assignment of the experimental time scales. 



V. CONCLUSIONS 

In this work, we have characterized the EET in the chloro- 
some by combining MD simulations, excited- state calcula- 
tions with TDDFT, and stochastic propagation of excitonic 
wave function. Our study shows that the coherent energy 
transfer occurs for 50 fs after the initial excitation, and the 
wavelike motion of the exciton is damped to vanish in 100 fs 
due to the dynamic disorder. The incoherent energy transfer 
subsequently takes place in the characteristic time constants 
from 1 ps to several tens of ps. The existence of multiple time 
scales reflects a hierarchy of the structures in the chlorosomes. 
We have interpreted the experimental time scales in terms of 
the theoretical results and also presented a possible verifica- 
tion of our assignment. 
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TABLE I: Results from MD/TDDFT calculations. Averages and 
standard deviations of the site energies, relaxation times of the site 
energies, the pure dephasing rates, and reorganization energies are 
provided for syn and anti monomers. 



(6)(cm-0 ^J{6^)(cm-^) TJfs) y(cm-^) ERJcm-^) 
syn 15243.0 532.2 8.88 938.8 658.7(266.5)^ 

anti 15319.4 556.5 4.51 525.8 719.9 (298.5^ 

^Reorganization energies calculated from t he sp ectral densities with 
the standard prefactor. [1331] 



TABLE II: Energy transfer time scales obtained from our simulations 
in comparison with the experimental time scales. 

Theoretical time scales 

Decoherence time 40-50 fs 

Equilibration over single roll 8 . 1 - 1 4. 1 ps 

Inter-layer transfer 1 . 1 - 1 .4 ps 

Inter-roll transfer 69.8 ps 

Fluorescence anisotropy decay 30 fs 
Experimental time scales 

One-color TA anisotropy decay (720-790 nm) [16] 1.7-3.9 ps 

One-color TA (720-790 nm) [.17] 200-300 fs, 1.7-1.8 ps, 5.4-5.9 ps, 30-40 ps 

Two-pulse photon echo (740-780 nm) \m} 140 fs, 1 ps, 10 ps, 280 ps 

Two-color TA (685 nm pump, 758 nm probe) [22] 120 fs, 1.1-1.2 ps, 12-14 ps, 46-52 ps 

2-D spectral line- shape dynamics [231 40 fs" 
" Decay time of the amplitude at the position of initial 
maximum f23ll . 



